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Abstract 



Oh 

-S 

In this article, we describe an approach for solving partial differ- 
ential equations with general boundary conditions imposed on arbi- 
trarily shaped boundaries. A function that has a prescribed value on 
the domain in which a differential equation is valid and smoothly but 
rapidly varying values on the boundary where boundary conditions 
qq are imposed is used to modify the original differential equations. The 

QO mathematical derivations are straight forward, and generically applica- 

ble to a wide variety of partial differential equations. To demonstrate 
the general applicability of the approach, we provide four examples: 
£\| (1) the diffusion equation with both Neumann and Dirichlet bound- 

t-H ary conditions, (2) the diffusion equation with surface diffusion, (3) 

the mechanical equilibrium equation, and (4) the equation for phase 
T". transformation with additional boundaries. The solutions for a few of 

these cases are validated against corresponding analytical and semi- 



analytical solutions. The potential of the approach is demonstrated 
with five applications: surface-reaction diffusion kinetics with a corn- 
er plex geometry, Kirkendall-effect-induced deformation, thermal stress 

in a complex geometry, phase transformations affected by substrate 
surfaces, and a self-propelling droplet. 



1 Introduction 

The smoothed boundary method [H EJ [3] and other similar approaches 
jH 13 16] have recently been demonstrated as powerful tools for solving various 
partial differential equations with boundary conditions imposed within the 
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computational domain. The method's origin can be traced to the embedded 
boundary method and the immersed boundary method (for an overview, see 
Ref . |3 [H El Uni E] ) • This method has been successfully employed in simu- 
lating diffusion processes [12] [13] and wave propagation [1] [2j [3j H2 US] con- 
strained within geometries described by a continuously transitioning domain 
indicator function (hereafter, the domain parameter) with a no-flux bound- 
ary condition imposed on the diffuse interface (as defined by the narrow 
transitioning region of the domain parameter). While those works demon- 
strated the potential for this type of numerical methods that circumvents 
the difficulties with constructing the finite element mesh (e.g., meshing the 
surface and then building a volumetric mesh based on the surface mesh or 
by combining regular subdomains that can be easily meshed), which is par- 
ticularly useful when dealing with complex structures. However, the method 
was only applicable to no-flux boundary conditions, and no approaches to 
extend the method to other types of equations or boundary conditions were 
available. Recently, a different formulation, based on asymptotic analy- 
ses, to solve partial differential equations in a similar manner was proposed 
[TBI [T7] [T8| [T9] 1201 I2T] . providing a justification of the method as well as 
increasing the applicability of the approach. 

In this paper, we provide a mathematically consistent smoothed bound- 
ary method and provide a precise derivation for the equations. The specific 
equations that we consider are: (1) the diffusion equation with Neumann 
and/or Dirichlet boundary conditions, (2) the bulk diffusion equation cou- 
pled with surface diffusion, (3) the mechanical equilibrium equation for lin- 
ear elasticity, and (4) Allen-Calm or Cahn-Hilliard equations with contact 
angles as boundary conditions. The method is especially useful for three- 
dimensional image-based simulations. 

2 Background 

The method is based on a diffuse interface description of different phases, 
similar to the continuously transitioning order parameters in the phase-field 
method [22] [23] [23] [25] [2BJ [27] often used in studying phase transformations 
and microstructural evolution in materials. In phase-field models, phases 
(which could be liquid, solid, vapor, or two different solids/liquids having 
different compositions) are described by one or more order parameters hav- 
ing a prescribed bulk values within each phase. In the interface, the order 
parameter changes in a controlled manner. Asymptotic analyses [27] can be 
used to show that the phase-field governing equations approach the corre- 
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sponding sharp interface problems in the sharp interface limit. 

We adopt this concept to describe internal domain boundaries by an 
order-parameter-like domain parameter, which may or may not be stationary 
and takes a value of 1 inside the domain of interest and outside. The 
equations will be solved where the domain parameter is 1, with boundary 
conditions imposed where the domain parameter is at the intermediate value 
(approximately 0.5). Figure [T] illustrates a schematic diagram of the sharp 
and diffuse interfaces. In the conventional sharp interface description, the 
domain of interest is O and is bound by a zero-thickness boundary denoted 
by 9f2 [Fig. [ija)]. Within f2, the partial differential equations need to be 
solved according to the boundary conditions imposed at dil. In the diffuse 
interface description, we employ a continuous domain parameter, which is 
uniformly 1 within the domain of interest and uniformly outside. In this 
case, the originally sharp domain boundary is smeared to a diffuse interface 
with a finite thickness indicated by < ij) < 1. Our target is to solve 
partial differential equations within the region where ip = 1 while imposing 
boundary conditions at the narrow transitioning interface region where < 
ij) < 1. By using this description, there is no specifically defined domain 
boundary. The system will determine the boundary by a variation of the 
domain parameter. In addition, the gradient of the domain parameter VV> 
will automatically determine the inward normal vector of the contour level 
sets of ij) (see Fig. [IJc)). 

3 Formulation 

3.1 General Approach 

The general approach is as follows. The domain parameter describes the 
domain of interest (-0 = 1 inside the domain, and ip = outside). The tran- 
sition between the two values described is smooth and taken as the solution 
to an Allen-Calm type dynamic equation (having a form of a hyperbolic 
function) described later. To derive the smoothed boundary formulation for 
Neumann boundary condition, the differential equation of interest (H) is 
multiplied by the domain parameter, ip. By using identities of the product 
rule of differentiation such as 



we obtain terms proportional to Vip. Since the unit (inward) normal of the 
boundary, ft, is given by V^/|Vi/>|, such terms can be written in dH/dn = 



ij)V 2 H = V • (ipVH) - Vij) ■ VH, 
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VH-n = V-ff-VVVI , and thus reformulated to be the Neumann boundary 
condition imposed on the diffuse interface. 

Similarly, to derive the smoothed boundary formulation for the Dirichlet 
boundary condition, the equation of interest is multiplied by the square 
of the domain parameter. Again using mathematical identities, ip 2 V 2 H = 
ipV ■ (ipVH) - i/jViJj ■ VH where ipVifj ■ VH = Vij) • V (i/)H) - H \Vip\ 2 , we 
obtain 

^ 2 V 2 H = ipV • (ipVH) - [Vip • V (ipH) - H\ VV'I 2 ] . (2) 

Note that H = H\qq associated with (V^l 2 appearing in the last term is the 
boundary value imposed on the diffuse interface. 

Specific details of the derivation depend on the equation to which the 
approach is applied, and we therefore provide four examples below. 



3.2 Diffusion Equation 

The first example is the diffusion equation with Neumann and/or Dirichlet 
boundary conditions. The Neumann boundary condition is appropriate, for 
example, as the no- flux boundary condition, while the Dirichlet boundary 
condition is necessary when the diffusion equation is solved with a fixed 
concentration on the boundaries. For Fick's Second Law of diffusion, the 
original governing equation is expressed as 

dC 

— = -V.j + S = V-(DVC) + S, (3) 

where j is the flux vector, D is the diffusion coefficient, C is the concen- 
tration, S is the source term, and t is time. Instead of directly solving the 
diffusion equation, we multiply both sides of Eq. ([3| by the domain param- 
eter if) that describes the domain of the solid phase: 

^^ t =^V-(DVC) + ^S. (4) 

Using the identity • {DVC) = V • (^jDVC) - Vip • (£>VC), Eq. Q 
becomes 

dC 

ip— = V • {tliDVC) - Vip ■ (DVC) + ipS. (5) 

Now, let us consider the boundary condition in this formulation. The Neu- 
mann boundary condition is the inward flux across the domain boundary, 
mathematically the normal gradient of C at the diffuse interface, and is 
treated as 

^ vv vy> • (DVC) dC 

n ' J = W\' J = W^T~ = ~ D aR = ~ Bm (6) 
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where n = V^/|V^| is the unit inward normal vector at the boundaries 
defined by the diffuse interface description. Equation ([6| can be rearranged 
to be Vip • (DVC) = Bn and substituted back into Eq. thus, we 
obtain 

dC 



V- 



dt 



V • (ipDVC) - B N + i/)S. 



To demonstrate that this smoothed boundary diffusion equation satisfies 
the assigned Neumann boundary condition (or specifying the boundary flux 
or normal gradient), we use the one-dimensional version of Eq. Q without 
loss of generality. By reorganizing terms and integrating over the interfacial 
region, we obtain 



if) 

a,i-£/2 



dt 



S\dx= 4>D 



dC 
dx 



a.i+Z/2 ra.i+i/2 
Oi-e/2 Ja.i-Z/2 



dip 
dx 



Bjydx, 



(8) 



where — £/2 < x < ai + ^/2 is the region of the interface, and £ is the 
thickness of the interface. Following Refs. [2j EJ [12j [15], we shall introduce 
the mean value theorem of integrals, which states that, for a continuous 
function, g(x), there exists a constant value, ho, such that: 



min g(x) < 



1 



q-p 



g(x)dx = ho < maxg(x), 



(9) 



where p < x < q. By eliminating the second term on the right-hand side 
of Eqs. ([7]) and Q, the no-flux boundary condition can be imposed; the 
resulting equation is similar to those proposed in Refs. [2j |3l [121 US]- How- 
ever, we retain the term in order to maintain the generality of the method. 
Therefore, the analysis presented here leads to an extension of the original 
method that greatly expands its applicability. 

Since the function on the left-hand side of Eq. ^ is continuous and finite 
within the interfacial region, we can use the mean value theorem of integrals 
to obtain the relation: 



a.i-Z/2 



dt 



S J dx = /io£. 



Using the conditions that ij) = 1 at x = ai + £/2 and ip = at x = ai 
the first term in the right-hand side of Eq. (|8|) is written as 



1 • D 



.dC 
dx 



0-D 



.dC 
dx 



D 



a, -5/2 



.dC_ 
dx 



Oi+C/2 



(10) 



(ii) 
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Since \Bip/dx\ = for x < a» — £/2 or x > en + £/2, the second term on the 
right-hand side of Eq. ([8]) can be replaced by 



9V 



9x 



Bjvdx 



<9V 



<9x 



Bj^dx. 



Substituting Eqs. (10), (11) and (12) back into Eq. QSJ) , we obtain 

"+0O 



ho£=D 



BC 
dx 



dtp 



dx 



Taking the limit of Eq. (13) for £ — ► 0, we obtain 

-+oo 



(12) 



(13) 



D 



BC 
dx 



5(x — a^Bjydx = B 



N 



(14) 

where BC '/8x\ ai+ £/ 2 = dC/dx\ ai and lim^_».o \dip/8x\ = 5(x — ai) when ip 
takes the form of a hyperbolic tangent function, and 8(x — Oj) is the Dirac 
delta function. The Dirac delta function has the property that J_°° S(x — 
a,i)f(x)dx = f(a,i), providing the second equality in Eq. (14). Therefore, 



Eq. ([14]) clearly shows that the smoothed boundary method recovers the 
Neumann boundary condition at the boundary when the thickness of the 
diffuse boundary approaches zero. This convergence is satisfied for both 
stationary and moving boundaries [12] . 

For imposing the Dirichlet boundary condition, we can manipulate the 
original governing equation in a similar procedure to the derivation of Eq. ([7| . 
Multiplying both sides of Eq. (15]) with ijj, we obtain 



,dc 

~3t 



V>V • (^DVC) - tpVtp ■ (DVC) + ip 2 S, 



(15) 



where the second term on the right-hand side can be replaced by ifi'Vip ■ 
(DVC) = D[Vip-V (V'C)-CV^-V^] = D[W-V(^C)-C|V^| 2 ]. Equation 
(|l~5l) is then rewritten as 



BC 

%b 2 — = i)V ■ {ipDVC) - D[Vyj • V (ipC) - C\Vip\ 2 } + ip 2 S, 



(16) 



where C in the third term will be the Dirichlet boundary condition, Bp, im- 
posed at the diffuse interface. Therefore, the smoothed boundary formulated 
diffusion equation with the Dirichlet boundary condition is 

BC 

t/; 2 — = VA7 ■ (i/jDVC) - DlVip • V (ibC) - B D \V^\ 2 ) + ^ 2 S. (17) 
at 
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To prove the convergence of the solution at the boundaries to the speci- 
fied boundary value, we again use a one-dimensional version of the smoothed 



boundary formulated equation. Integrating Eq. (17) over the interfacial re- 
gion and reorganizing terms, we obtain 



ai+4/2 



,dC 
~dt 



d_ 

dx 



ipD 



OX J J Jai-i/2 \ dx J L dx 

(18) 

Similar to the derivation of Eq. (10), the left-hand side of Eq. (18) is pro- 



B D 



dip 
Ox 



dx. 



portional to the interfacial thickness and approaches zero in the limit of 
£ —* 0. On the right-hand side of Eq. ( 18 ), the gradient of ip approaches the 
Dirac delta function, 8{x — ai), as the interface thickness approaches zero. 



Therefore, we can reduce Eq. (18) to 



= 23 



dipC g dip 
dx dx 



dipC 
dx 



B D 



dip 
dx 



(19) 



in the limit £ — > 0. By integrating over the interfacial region of Eq. (19) 
again, we obtain 



l-C 



o-c 



which in the limit of £ — > gives 



aiH/2 Qlp , 
B^—dx, 

oj-e/2 



ai+^/2 



+oo 



(5(x — a,i)Br>dx = 29^ 



(20) 



(21) 



Therefore, the smoothed boundary formulation recovers the specified Dirich- 
let boundary condition: C\ ai = Bjj\ ai . 

In this method, the boundary gradient, Bn, and the boundary value, 
Bjj, are not specified to be constant values. They can vary spatially and/or 
temporally or be functions of C or other parameters. In addition, one can 
impose Neumann and Dirichlet boundary conditions simultaneously to yield 
mixed (or Robin) boundary conditions. The equation then becomes 



iP'- 



,dC 

Ik 



ipV-(ipDVC)-ip\Vip\ N B N (x)-V^-D[V(ipC)-B D (x)ViP] D +iP 2 S, 

(22) 

where Bn(x) and Bd(x) are spatially dependent Neumann and Dirichlet 
boundary conditions specified at different parts of the boundary, and the 
subscripts W and '23' denote the quantities associated with the boundaries 
to which the Neumann and Dirichlet boundary conditions are imposed. 
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3.3 Surface Diffusion Coupled Bulk Diffusion 

The second example will demonstrate that surface diffusion can be imple- 
mented into the smoothed boundary equation derived above. For this case, 
we take the set of equations that includes surface reaction, bulk diffusion 
and surface diffusion to describe an oxygen reduction model in a solid oxide 
fuel cell (SOFC) cathode [28J. The oxygen vacancy concentration, C, on the 
cathode surface is governed by Fick's Second Law: 

where n, s and r are the unit normal, primary tangent and secondary tangent 
vectors of the surface, respectively. Here, the parameter I is the characteris- 
tic thickness of the surface and is multiplied into the surface Laplacian term 
to maintain the dimensional agreement of the equation. The parameters 
D b , At, D s , L and t are the bulk diffusivity, reaction rate, surface diffusivity, 
accumulation coefficient and time, respectively. Thus, the term on the left- 
hand side represents the flux from the bulk, and the terms on the right-hand 
side represent the surface reaction, surface Laplacian and concentration ac- 
cumulation, respectively. For simplicity, these parameters are all assumed to 
be constant. In the bulk of cathode particles, the oxygen vacancy diffusion 
is also governed by Fick's Second Law: 

^ = D b V 2 C. (24) 
To simulate the oxygen vacancy concentration evolution in the cathode, 



the two diffusion equations, Eqs. (23) and (24), are coupled and need to 



be solved simultaneously. In this case, the two equations will share the flux 
normal to the cathode surface as the common boundary condition. Recently, 
this set of equations was formulated using the concept of diffuse interface 
approach [19] , which leads to two differential equations that are coupled by 
boundary conditions. We will show below that the coupling can be achieved 
by applying the smooth boundary formulation described herein to obtain 
one single equation that governs both surface and bulk effects. 



The derivation is as follows. We first multiply Eq. (24) with ip and 
applying the product rule of differentiation to obtain 

BC 

^~dt = ° bV ' ^ VC) ~ DbVi) ' VC ' (25) 
As in Eq. (|6|, the normal derivative to the diffuse interface is defined by 



dC/dn = VC • Vip/\Vip\. Substituting this relation back into Eq. (23) and 
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rearranging terms give 



MA 

D h 



kC - ID, 



a 2 d 2 \ T dc 







Ot 



Substituting Eq. (26) into the second term in Eq. (25), we obtain 

d 2 



dC 

^— = D b V ■ (^VC) - |VV| 



kC - ID, 



d 2 \ c 

ds 2 dr 2 J 



dC_ 
~dt 



(26) 



(27) 



This equation combines the bulk diffusion and surface diffusion into one sin- 



the interfacial thickness approaches zero, Eq. (27) will reduce to Eq. (23) at 



gle equation, and will be used in examples presented in Sections 4.2 and 5.1 
In the bulk (|W>| = and ip = 1), Eq. d27} reduces back to Eq. pih. When 



the interface (iVf/'l 7^ 0) as has been proven in Section 3.2 



To calculate the surface Laplacian, we use the following method. The 
unit vector of the concentration gradient is given by p = VC/|VC|. The 
unit secondary tangential vector on the surface can be obtained by f = 
(n x p)/\n x p\, and the unit primary tangential vector is then obtained 
by s = (r X n)/\T x n\. In this case, the surface flux has no projection 
in the r direction {p ' ■ f = 0). We can calculate the surface diffusion flux 
along the primary tangent direction simply by projecting the concentration 
gradient into the primary tangential direction. The surface flux is calculated 
by taking the inner product between the concentration gradient and the unit 
primary tangential vector for magnitude, and it is along the opposite primary 
tangential direction: 

1 = -W s (p- s)s. (28) 

Since the Laplacian operator is independent of the selection of coordinate 
system, the value of the surface Laplacian can be then obtained by taking 
the negative divergence of the surface flux: 



W, 



ds 2 dr 2 



C 



-V • j s 



(29) 



where V • j s is the divergence of j s in the global Cartesian grid system of 
the computational box. 

To simulate only the surface diffusion on a diffuse-interface described 
geometry, one can simply eliminate all bulk-related terms to obtain 

dC 



L 



dt 



-V • j s 



(30) 



such that only a concentration evolution along the interfacial region will 
occur. 
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3.4 Mechanical Equilibrium Equation 



The smoothed boundary method can also be applied to the mechanical equi- 
librium equation. When a solid body is in mechanical equilibrium, all the 
forces are balanced in all directions, as represented by 



do 



0, 



(31) 



where the subscript 'i' indicates the component along the i-th direction, and 
<Jij is the stress tensor. Repeated indices indicate summation over the index. 
For a linear elasticity problem, the stress tensor is given by the generalized 
Hooke's Law: 

&ij = Cijki(eki — pSki), (32) 

where is the elastic constant tensor, and p is a scalar body force, such 
as thermal expansion (aAT) or a misfit eigen-strain (e° = (a p — a m )/a m ), 
which depends on the governing physics. The total strain tensor is defined 
by the gradients of displacements as 



1 / dui duj 

2 \dxj dxi 



(33) 



where u% is the infinitesimal displacement along the i-th direction. Substi- 



tuting Eqs. (33) and (32) back into Eq. (31) gives 
d 



c 1 ( Out du^ 

dxj 8,3 2 I dxi dxk 



d 

d Xj 



pCijkiSki 



(34) 



We can multiply Eq. ( |34[ ) by the domain parameter that distinguishes the 
elastic solid region (■0 = 1) from the environment (ip = 0) to perform the 
smoothed boundary formulation. After collecting the terms associated with 
dij)/dxj on one side of the equation, we obtain 



_d_ 

dxj 



ijkl 



1 f duk dui ^ 



2 \ dxi dxk J 



dip 
dxj 



ijkl 



1 f duk chM 

2 V dxi dx k 

- JL 
dxi 



pS, 



id 



ippCijkiSki 



(35) 



The traction exerted on the solid surface is defined by iVj = —a^rij, 
where rij is the inward unit normal of the solid surface. We again use the 
definition of the inward unit normal of the boundary: ni = Vt/>/|V?/>|- (In 
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the indicial notation, d'p/dxi = Vip and \J (dip /dxi)(dip /dxi) 
Therefore, the traction force is given by 



Nr 



ijkl 



du k du^ 
dxi dx k 



PS. 



|v#) 



(36) 



Substituting Eq. (36) back into Eq. (35) returns the smoothed boundary 
formulation of the mechanical equilibrium equation with a traction boundary 
condition on the solid surface: 



_d_ 

dxj 



1>G 



ijkl 



1 / duk Qui \ 

2 1 dxi dx k J 



+ \Vi/>\Ni 



_d_ 

dxj 



ippCijkiSki 



(37) 



where d(tppCijki5ki)/dxj = pi can be treated as an effective body force along 
the i-th direction. 

For linear elasticity problems with presribed displacements at the solid 
surface, one can perform the smoothed boundary formulation as in the 
derivation of the Dirichlet boundary condition in Section [3.2| by multiplying 
Eq. (34) by ip 2 and using the product rule to obtain 



d 

dxj 



1 



+ 



dm 



( duk 
2 \ dxi ' dx 



dip 
dxn 



a 



dip\ 
dxj 



Cii 



ijkl ; 



dip dip 
dxi dx k 



ijkl 



r 



( dipu k dipuA 



1 

2 V dxi 
d 



dx k J 



dxi 



pCijkiS, 



ijklVkl I j 



(38) 



where the displacements Uk and ui appearing in the third term on the left- 
hand side will be the boundary values of the displacements at the solid 
surface. An equivalent formulation for the mechanical equilibrium equation 
can also be obtained by asymptotic approach 



3.5 Equations for Phase Transformations with Additional 
Boundaries 

Phase transformations affected by a mobile or immobile surface or other 
boundary are of importance in many materials processes including hetero- 
geneous nucleation that takes place at material interfaces |3U|I31|. Maintain- 
ing a proper contact angle at the three-phase boundary (where the interface 
between the two phases meets the surface) is necessary in capturing the 
dynamics accurately, since the contact angle represents the difference be- 
tween surface energies (tensions) of different phase boundaries. While there 
are previous works that developed a method to impose the contact-angle 
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boundary condition |30| I3T] on sharp domain walls, here we show that a 
similar model with diffuse domain walls can be obtained simply by apply- 
ing the approach described above. Below, we assume that the boundary 
is immobile, but this assumption can be easily removed by describing the 
evolution of the domain parameter as dictated by the physics of the system. 

In the Allen-Cahn and Cahn-Hilliard equations of the phase field model, 
the total free energy has the following form [221 E3] : 



F 



dfi, (39) 



where <f> is referred to as the phase field or order parameter commonly used 
to define different phases, and e is the gradient energy coefficient in the 
phase field model. We take the variational derivative according to Euler's 
equation: 



SF = J f||-e 2 V 2 ^dO + J (e 2 V(j>-fi)dA, 



(40) 



where ft is the unit normal vector to the domain boundary d£l. The bulk 
chemical free energy, /, is a double-well function of eft. (This can also be 
derived from the order parameter <fi changing with local "velocity" <f>.) For an 
extremum of the functional F, SF = must be satisfied. This requirement 
provides two conditions: 

|^ - e 2 V 2 = in O, (41a) 



e 2 X7(f> ■ ft = on dn. 



(41b) 



Following Eq. (41a), we find 



W = eVW = -V(|W), 



(42) 



which can be rewritten as V/ = V(e 2 | V(j)\ 2 )/2. We thus find a useful equal- 
ity for deriving the contact angle boundary condition: 



2/ 



(43) 



In the smoothed boundary method, we introduce a domain parameter tjj 
to incorporate boundary conditions in the original governing equation. As 
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mentioned earlier, the level sets of this domain parameter tp describe the 
original boundaries and should satisfy n = Vip/\Vip\. On d£l, we impose a 
contact angle 9. Thus, 



n 



Vtp Vcf) 



|V0| |V^| |V0| 



cos 9. 



(44) 



Substituting Eq. (43) into Eq. (44), one derives the following boundary 

/2f 



condition formulation: 



Vip-V<p = |V^| cos9 



(45) 



This contact-angle boundary condition is similar to the one suggested by 
Warren et al. [31J for contacting a sharp interface, for which a Dirac delta 
function will replace |V^|. 

The bulk chemical potential is defined by the variational derivative of 
the total free energy of the system: 



li 



_ df_ 

Sep dcf) 



e 2 V 2 c 



(46) 



as it appeared in the first term of Eq. ( 40 ) . Multiplying both sides of Eq. ( 46 ) 
by the domain parameter tp gives 



ipfj, = tp 



df 



iPe 2 V 2 (t> = ip 



df 



e 2 V • (iN<t>) + e 2 Vtp ■ Vcp. (47) 



We substitute the contact- angle boundary condition, Eq. (45), into the third 



term in Eq. (47) and obtain the smoothed boundary formulation for the 



chemical potential by dividing both sides by ip: 
df e 2 , e|W| 



2/cos6». 



(48) 



For a nonconserved order parameter in the phase field models, the evo- 
lution is governed by the Allen-Cahn equation (2¥], in which the order pa- 
rameter evolves according to the local chemical potential variation: 



at 



-Mfi 



(df_ _ 
\d(f) ip 



V • (^V0) + 



2/cos<9 



(49) 



For a conserved order parameter, the evolution of the order parameter is 
governed by the divergence of the order-parameter flux, while the flux is 
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proportional to the gradient of the chemical potential gradient. This process 
is governed by the Cahn-Hilliard equation |22| [25] : 



^ = V-(MV M ), (50) 
for which the smoothed boundary formulation is obtained by (see Section 



3.2) 



ip-^ = V • (ipMVfi) - Vip • (MV/i). (51) 

Note that — MVjU = j is the flux of the conserved field order parameter. 
Therefore, the second term represents the fluxes normal to the domain wall 
(equivalent to Eq. Q): 

V%b-{MVfi) = -(j-n)\V^j\. (52) 

Substituting the flux across the domain wall, our final smoothed boundary 
formulation of the Cahn-Hilliard equation is then written as 



ot 



^v(|-^V.(^) + ^v^7cos. 



+ \V4>\J n , (53) 



where J n = j ■ n. In practice, ip has a very small cutoff value such that the 
terms containing can be numerically evaluated. For time dependent 
problems, the equation is divided by ip before numerical implementation. 



4 Validation of the approach 

We demonstrate the validity and accuracy of the approach using bulk/surface 
diffusion in Sections |3.2|and 3.3 as well as the phase transformation of three 



phase systems in Section 3.5 



4.1 ID Diffusion Equation 

First, we perform a ID simulation to demonstrate that the Neumann and 
Dirichlet boundary conditions are satisfied on two different sides of the do- 
main. Fick's second diffusion equation, Eq. ([3]), with the given source term 
is solved within the solid phase that is defined by ip = 1. The diffusion coef- 
ficient D is set to be 1, and the source S is 0.02. On the right boundary of 
the solid, the gradient of C is set to be -0.05, while on the left boundary, the 
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value of C is set to be 0.4. We perform the smoothed boundary formulation, 
as in the derivation of Eq. (22), to obtain 



V' 2 ^ = V'V-^VC)-^[|VV|(-0.05)] r -[VV-V(VC)-|VV| 2 (0.4)] i +^ 2 (0.02), 

(54) 

where the subscripts V and T indicate the right and left interfaces. The 
solid region is approximately in the range between the 102-th and 298-th 
grid points. We use a hyperbolic tangent form for the continuous domain 
parameter ifi as 

ip = ^{tanh [0.8(x - 10) + 1] - tanh [0.8(x - 30) + 1]}, (55) 

such that the interfacial thickness is taken to be approximately 6 grid spac- 
ings. The initial concentration is C = everywhere in the computational 
box. A standard finite central differencing scheme in space and the Euler 
explicit time scheme are employed in the simulation. The grid spacing is 
taken to be Ax = 0.1. 

Figure [2] shows the concentration profiles taken at four different times 
(in blue solid lines). The domain parameter is plotted in the red dashed 
line. On the right interface, it can be clearly observed that dC/dx = —0.05 
at all times, except for a rapid change from dC/dx = to dC/dx = —0.05 in 
the very early transient period. In the early period, the concentration even 
takes negative values to satisfy the gradient boundary condition imposed at 
the right boundary. On the other hand, the concentration remains at 0.4 
during the entire diffusion process, except in the very early transient period 
during which C changed from to 0.4. This result clearly demonstrates 
both Neumann and Dirichlet boundary conditions are satisfied on the diffuse 
interfaces. 



4.2 Surface Diffusion and Bulk Diffusion in a Cylinder 

To further demonstrate the validity of the smoothed boundary method, we 
apply the method to a cylinder for which a cylindrical coordinate grid system 
can be used. We solve the coupled surface-bulk diffusion problem using 
both the smoothed boundary and standard sharp interface formulations in 
the same grid system for comparison. Again, we use a continuous domain 
parameter ip to define the solid region of a cylinder (ip = 1 for solid, and 
ip = for environment). The solid surface is then represented by < ip < 1. 



For the smoothed boundary case, we solve Eq. (27) using the central 



differencing scheme in space. The radial direction is discretized into 80 grid 
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points, and the longitudinal direction is discretized into 600 grid points. 
The grid spacing Ax is 1/60, such that the radius of the cylinder is R = 
60Ax = 1, and the length of the cylinder is 10i?. The thickness of the 
diffuse interface is approximately 4~5 grid spacings; thus, the characteristic 
thickness appearing in Eq. (27) is set to be Z = 4.5 Ax = 0.075. Here, we set 
the surface accumulation coefficient L to be for simplicity. We investigate 
two cases: one with a low surface reaction rate, k = 2.1, and the other with 
a high surface reaction rate, k = 1000. 

To compare the results, we solve the original form of the coupled sur- 
face and bulk diffusion equations using the sharp interface approach with 
the same finite difference method. The same grid system is used, except a 
cylinder surface is now explicitly placed at R = 1 at which the boundary 
condition is imposed. In this case, we calculate the normal flux to the sur- 
face by the right-hand side of Eq. (23) with the grid system on the cylinder 
surface, and then use the flux as the boundary condition for the concen- 
tration evolution in the bulk and on the surface, Eq. (24). Note that the 
characteristic thickness / drops in the sharp interface description as the limit 
/ — > is taken. 

Figures [3^a) and (b) show the concentration profiles in the cylinder at 
the steady state obtained using the smoothed boundary method and the 
finite central difference method. For clarity, only the concentration in the 
region of < z < 5R is presented. The top rows in Figs, ga) and (b) arc 
the smoothed boundary results, and the bottom rows are the sharp interface 
results. The results from the two methods are clearly in excellent agreement. 

Shown in Figs. |4ja) and (b) are the concentration profiles plotted along 
the longitudinal line at r = 0, r = 2R/3, and r = R, respectively. Again the 
plots show that the differences between the results from the two methods 
are small for the cylindrical geometry. As mentioned in a previous section, 
the error of the smoothed boundary method is proportional to the interfa- 
cial thickness. Based on our tests, we found that, even for an interfacial- 
thickness-to-radius ratio of around 1/5, the maximum error between the two 
methods appearing near the surface is still around 2% (shown in the solid 
square markers), while the error in the bulk region is significantly smaller 
than that. If we select the interfacial-thickness-to-radius ratio to be 1/10, 
the maximum error appearing in the entire solid region is on the order of 
1 x 10~ 3 (including the region near the surface). Another controlling factor 
of errors is the number of grid points across the diffuse interface. From our 
numerical tests, we noticed that at least 4 grid points are required to prop- 
erly resolve the sharp change in i/j across the interface such that the errors 
are reasonably small. In addition to the steady state solution, the transient 



16 



state solutions are also in excellent agreement during the entire diffusion 
process. This demonstrates that the smoothed boundary method can be 
employed to accurately solve coupled surface diffusion and bulk diffusion 
problems. 



4.3 Contact Angle Boundary Condition 

We perform a simple 2D simulation to validate the smoothed boundary 
formulation for the contact-angle boundary condition at the three-phase 



boundary. Equations (49) and (50) are tested for nonconserved and con- 
served field order parameters, respectively. The computational box sizes are 
L x = 100 and L y = 100, and the parameters used are Ax = 1, M = 1, and 
e = 1. On the computational box boundaries, the normal gradients of the 
order parameter are set to be zero: d<p/dx = at x = and x = 100 and 
d(f)/dy = at y = and y = 100. A horizontal fiat wall is defined by a 
hyperbolic tangent function of the domain parameter tp 

ij> = 1 tanh (y - 30) + ^, (56) 

such that tfj = 0.5 is at y = 30 and tp gradually transitions from to 1 
from below the wall to above. The wall thickness is approximately 5 grid 
spacings. The initial phase boundary is vertically placed at the middle of 
the domain (x = 50) with phase 1 (if) = 1) and phase 0(^ = 0) on the left 
and right halves, respectively. 



In the first case with nonconserved order parameter, we evolve Eq. (49) 
with a 60-degree contact angle. The result clearly shows a 60-degree contact 
angle at the three-phase boundary as imposed (Fig. [5^a)). The angle can 
be measured by the intersection between the two contours of ip = 0.5 and 
4> = 0.5, as shown in Fig.[5^b). The 60-degree angle is maintained during the 
entire evolution, except for the very early transient period when the contact 
angle changed from 90 to 60 degrees. Due to the imposed contact angle, the 
initially fiat phase boundary bends and creates a negative curvature of phase 
1. As a result, the phase boundary moves toward phase 0, and eventually 
only phase 1 remains in the system. 



For the second case with conserved order parameter, we evolve Eq. (50) 
with a contact angle of 120 degrees. As expected, the phase boundary inter- 
sects the wall at a 120-degree contact angle (Fig. |5jc) and (d)). In contrast 
to the Allen-Cahn type dynamics, due to the conservation of the order pa- 
rameter, the phase boundary near the wall moves toward the left while the 
phase boundary away from the wall moves in the opposite direction. As a 
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result, the phase boundary deforms to a curved shape. When the system 
reaches its equilibrium state, the phase boundary forms a circular arc with 
a uniform curvature everywhere along the phase boundary, such that the 
total surface energy is minimized (see Fig. pfc) for t = 1.3 x 10 5 ). 



5 Applications 

While the details of the scientific calculations performed by applying of these 
methods will be published elsewhere, it is worthwhile to show some of the 
results to demonstrate the potential of the method. 



5.1 Surface- Reaction Diffusion Kinetics 

The first example is ionic transport through a complex microstructure. Here, 
the ion diffusion is driven by a sinusoidal voltage perturbation. For the 
steady state solution, the time dependence of the form exp(iu;i), where u 
is the angular frequency and i = y/— T, can be removed, as in the equation 
derived by Lu et al. [28 . For a demonstration, we solve the steady state 
solution for the case without surface diffusion while solving the transient 
state solution for the case with surface diffusion. For the first case, the 
smoothed boundary formulated equation is given by 

V- {ipVC) - \ViP\kC = iujipC, (57) 

where C is the concentration amplitude consisting of a real and imaginary 
part. This equation is solved by a standard alternative direction iterative 
(ADI) method in a second-order central-difference scheme in space (Ax = 
0.04). For the transient state solution, we keep the time dependence as 



is, and the smoothed boundary formulation is given by Eq. (27) in which 
surface diffusion, bulk diffusion and surface reactions are all considered. For 
simplicity, the surface accumulation term is ignored (L = 0). Equation 
( |27| ) is solved by a second-order central-difference scheme in space (Ax = 
0.04) and the Euler explicit time stepping scheme (At = 0.01). Here, we 
employed an Allen-Cahn type equation [32], |33l [M] to smooth the initially 
sharp boundaries of experimentally obtained 3D voxelated data (ip = 1 for 
the voxels in the cathode and ip = for the voxels in the pores): 

at = -ty + 6 v ^ " eV2f WW x ' (58) 

where / = ip 2 (l — ip) 2 is a typical double- well function, and e is the gradient 
energy coefficient. The interfacial thickness (0.1 < ip < 0.9) is given by 
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2ey/2. Note that the third term in Eq. (58) is used to remove the curvature 
effect such that the location of ip = 0.5 does not change during the smoothing 
process if x = 1- The computational box contains 321, 160 and 149 grid 
points along the x, y and z directions. 

Figure [6] shows the steady-state concentration for the case in which 
Df, = 1, k = 0.1, D s = and to = 0.55. The boundary condition on the com- 
putational box is C = 1 at y = 0, C = at y = 6.4 and no-gradient on the 
remaing four sides. As shown in Fig. |6^a), the real part of the concentration 
decays from 1 to over the complex cathode microstructure to satisfy the 
boundary condition given on the domain box at y = and y = 6.4. For the 
imaginary part, the values at y = and y = 6.4 remain at as assigned. In 
the middle region, a negative value of the imaginary part occurs due to the 
phase shift resulting from the delayed response. 

Figure [7] shows the concentration distribution taken at two different times 
for the case in which = 1, k = 2.1 and D s = 10, with DC loading 
(u = 0). The boundary conditions on the computational box are the same as 
in the AC loading case above. An enhanced concentration along the irregular 
surface due to surface diffusion can be clearly observed in the intermediate 
stage, Fig. [7](a). As the concentration propagates through the bulk region, 
the system eventually approaches its steady state, and the concentration 
enhancement diminishes, Fig. [T^b). Figures[7^c) and (d) are magnified views 
of (a) and (b). 

The smoothed boundary method can also be used to impose Dirichlet 
boundary conditions on irregular surfaces. For example, if the ion diffu- 
sivity is much higher in the electrolyte phase than in the cathode phase, 
the concentration in the electrolyte will be nearly uniform. To simulate 
this scenario, we impose a fixed concentration at the electrolyte-cathode 
contacting surface as the boundary condition. On the computational box 
boundaries, we set C = at y = 10.44 and the no-flux boundary condition 
for the remaining five sides. The material parameters are selected to be 
Df, = 1, k = 0, and D s = 0. Figure [8] shows the simulation results for a pure 
bulk diffusion example with a fixed value C = 1 imposed at the LSC (cath- 
ode) -YSZ (electrolyte) interfaces. In this case, since the contacting areas 
are small (compared to the cross-sectional area of LSC on the x-z plane in 
Fig. [6|a)), ion diffusion along the lateral directions (x and z) is large. As 
a result, the concentration drops very rapidly in a short distance from the 
contacting areas. Therefore, the concentration distribution is very different 
from the ones shown in Figs. [6] and [7J where the cross-sectional areas at 
y = and y = 6.4 (x-z planes on the computational box) are approximately 
equal. 
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5.2 Kirkendall Effect Diffusion with a Moving Boundary Driven 
by Coupled Navier-Stokes-Cahn-Hilliard equations 

The third application will demonstrate the smoothed boundary method's 
broad applicability by applying it to the coupled Navier-Stokes-Cahn-Hilliard 
equations [351 EU ETJ [381 [39] . This particular formulation aims to solve dif- 
fusion problems with the Kirkendall effect with vacancy sources and sinks 
in the bulk of the solid |40[ Rd| S2J 03] E]. In this case, the solid experiences 
deformation due to vacancy generation and elimination. The Navier-Stokes- 
Cahn-Hilliard equations are coupled to the smoothed boundary formulation 



of the diffusion equation in Section 3.2 as a model of plastic deformation 



due to volume expansion and contraction resulting from vacancy flow. 

When the diffusing species of a binary substitutional alloy have different 
mobilities, the diffusion fluxes of the two species are unbalanced, creating 
a net vacancy flux toward the fast diffuser side. Here, we denote the slow 
diffuser, fast diffuser and vacancy by A, B and V, respectively. Due to the 
accommodation/supply of excess/depleted vacancies, the solid locally ex- 
pands / shrinks |45[ |4"6] |4"T] 08] when maintaining the vacancy mole fraction 
at its thermal-equilibrium value. We treat the solid as a very viscous fluid 
\4Q\ [50] IST] [52] with a much larger viscosity than the surrounding environ- 
ment. In this case, we solve the Navier-Stokes-Cahn-Hilliard equations to 
update the shape of the material as follows [53] : 

,(2n r \ i 



-VP + V- (r/Vv) +V( -jSvj +— /xV^ = 0, (59a) 



V-v = -5y, (59b) 



^-W^MV 2 ^-^, (59c) 

where P is the effective pressure, r\ is the viscosity, v is the velocity vector, 
d is the number of dimension, and C a is the Cahn number reflecting the 
capillary force compared to the pressure gradient. One great convenience of 
solving this type of phase field equation is that it automatically maintains 
the domain parameter in the form of a hyperbolic tangent function while 
updating the location of the diffuse interface. Note that we have ignored 
the inertial force in the Navier-Stokes equation to obtain Eq. (59a) since 



the deformation is assumed to be a quasi-steady state process. The vacancy 
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generation rate that results in the local volume change is given by 

V • (D VB VX B ) 

<V = JZ Y eq, , (60) 

where Xb is the mole fraction of the fast diffuser, Xy is the thermal- 
equilibrium vacancy mole fraction, Dvb is the diffusivity for vacancy flux 
associated with VXg, and p\ is the lattice site density. The fast diffuser mole 
fraction evolution is governed by the advective Fick's diffusion equation as 

<jA '" vVXb = V-(D^ b VX b )-X b S v , (61) 



dt 



where D^ BB is the diffusivity for a fast diffuser flux associated with VXg, and 
the advective term accounts for the lattice shift due to volume change. Since 
the diffusing species cannot depart the solid, a no-flux boundary condition 
is imposed at the solid surface. Thus, we obtain the smoothed boundary 



formulation of Eq. (61) as 

^yW 1 ~ V ' VXb ) = V ' W D bbVXb) - ^X B S V . (62) 
As the concentration evolves, the shape of the solid is also updated by 



Eq. (59c) and iteratively solving Eqs. (59a) and (59b) by applying a projec- 



tion method jSH [55] . 

The slow and fast diffusers are initially placed in the left and right halves 
of the solid, respectively. We use theoretically calculated diffusivities for this 
simulation [561 E3 Ell [59] . Figure [9] shows snapshots of the concentration 
profiles (left column) and velocity fields (right column) from a 2D simulation. 
As the fast diffuser diffuses from the right to the left side, the vacancy 
elimination and generation cause contraction and expansion on the right and 
left sides, respectively. As a result, the initially rectangular slab deforms to 
a bottle-shaped object. 

In another scenario in which the vacancy diffusion length is comparable 
to or smaller than the distance between vacancy sources and sinks, the 
explicit vacancy diffusion process must be considered 59. 60 . In this case, 
vacancies diffuse in the same manner as the atomic species. In the bulk 
of a solid devoid of vacancy sources /sinks, the concentration evolutions are 
governed by 

£ ^± = V-(D VV VX V + D VB VX B ), (63a) 
at 

^ = V • {D BV VX V + D V BB VX B ). (63b) 
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Since the solid surfaces are very efficient vacancy sources / sinks |53| 161 j , we 
impose the thermal-equilibrium vacancy mole fraction at the solid surfaces 



as the Dirichlet boundary condition for solving Eq. (63). In this case, the 



smoothed boundary formulation of Eq. ( 63 ) is given by 



= ^ V ' ^( D wVX v + D VB VX B )} - K, (64a) 



^^W = ^ V ' M D *v VX v + D bbVX b )} + Y^epK, (64b) 



where K = Dyy\S/i[) • V(ipXy) — \Vip\ 2 Xy]. Since the vacancy generation 
and elimination in this scenario only occurs on the solid surfaces, no internal 
volume change needs to be considered in the bulk. Therefore, instead of 
using a plastic deformation model as in the previous case, we adopt a typical 
Cahn-Hilliard type dynamics to track the shape change: 



dt M |v^| i - x v 



where Jy = DyyVXy + Dys^Xg is the vacancy flux, and the last term 
represents the normal velocity of the solid surfaces due to vacancy injection 
into or ejection from the solid. 

An example of the results obtained using this approach is the growth 
of a void in a rod (62J |63l El]- The above equations are solved using a 
central difference scheme in space and an implicit time stepping scheme. The 
vacancy mole fraction is fixed at the void and cylinder surfaces. The fast 
diffuser initially occupies the center region, while the slow diffuser occupies 
the outer region. A void is initially placed off-center in the fast diffuser 



region. Figure 10 shows snapshots of the fast diffuser mole fraction profile 
and the vacancy mole fraction profile (normalized to its equilibrium value). 
As the fast diffuser diffuses outward, vacancies diffuse inward from the rod 
surface to the void surface, causing vacancy concentration enhancement and 
depletion in the center and outer regions, respectively. To maintain the 
equilibrium vacancy mole fraction at the rod and void surfaces, vacancies are 
injected and ejected at those surfaces. As a result, the rod radius increases, 
and the void grows. Such dynamics was examined using a sharp interface 
approach [61, but this new method provides the flexibility in geometry to 



examine cases where a void initially forms off-center. 
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5.3 Thermal Stress 



Since an SOFC operates at temperatures near 500~1000°C, the thermal 
stress is important for analyzing mechanical failure. Here, we expand the 



generalized mechanical equilibrium equation, Eq. (37), for a linear, elastic 



and isotropic solid. In this case, the elastic constant tensor is expressed by 

An = Cnii = C2222 = C3333) (66a) 



A 
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1133, 



(66b) 



A44 =Cl212 = Cl221 = C2H2 = C*2121 = C2323 = C2332 
=C3223 = C3232 = C1313 = C1331 = C3113 = C3131. 



(66c) 



The remainder of the elastic constant components vanish. We use the co- 
ordinate notation to replace the indices i = 1, 2 and 3 by x, y and z, 
respectively. The infinitesimal displacements along the x, y and z directions 



are then replaced by u, v and w, respectively. Thus, Eq. (37) in the three 
Cartesian directions is rewritten as 
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|V^|iV 2 
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dz 



[Vp(Aii + 2Ai 2 )], 



for the x, y and z directions, respectively. To numerically solve this equation, 



we reorganize the terms in Eqs. (67a)-(67c) to form 



C\u = h\, C2V = h,2, and C3W = /13 



(68) 



where £\, £2 and £2 are the linear differential operators associated with 
u, v and w, respectively, in Eq. (67). 
/13, are the remaining terms collected in Eq. (67). 



The right-hand sides, h%, h^ and 
The linear differential 



operators are discretized in the second-order central differencing scheme in 
space. We employ an ADI solver for the linear differential operators and 



iterate Eq. (68) until u, v and w all converge to their equilibrium values. 

We select the material parameters as follows: a YSZ AT = 1% and 
a LSC AT = 2%. The elastic constants are choosen arbitrarily as Xjf z = 
20 x 10 7 , \22 SZ = 10 x 10 7 , and \%f z = 5 x 10 7 (dimensionless) such 
that the solid is isotropic in mechanical behavior, (An — Ai2)/(2A44) = 1. 
The LSC (cathode) phase is softer than the YSZ (electrolyte) phase, and 
its elastic constant is assumed to be 0.75A^ . Again, we use domain pa- 
rameters to indicate the YSZ phase (tpysz = 1 inside YSZ and i/jysz = 
outside YSZ) and LSC phase (tpisc = 1 inside LSC and ipLSC = out- 
side LSC). The entire solid phase is then indicated by the sum of the two 
phases, ip = ipysz + tpLSC = 1- The body force term and elastic constant 



tensor in Eq. (67) are replaced by an interpolated, spatially dependent ther- 
mal expansion, tpp = O.Olipysz + 0.02ipLSC, and elastic constant tensor, 
Ajj = A« ipYSZ + ^ff C,l pLSC- The solid surface is assumed to be traction- 
free, Ni = 0. 

In this simulation, we select a computational box containing 160, 160, 
and 149 grid points along the x, y, and z directions, respectively. The grid 
spacing is Ax = 0.04. We assume a rigid computational box with frictionless 
boundaries on the six sides, which means that u = 0, v and w are free on the 
two y-z planes, v = 0, u and w are free on the two x-z planes, and w = 0, u 
and v are free on the two x-y planes of the computational box boundaries. 

Figure [TT|a) illustrates our experimentally obtained microstructure con- 
taining the cathode (LSC) and electrolyte (YSZ) phases. The yellow color 
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indicates the cathode phase, and the cyan color indicates the electrolyte 
phase. Shown in Fig. [XT|b) are the calculated Von-Mises stresses resulting 
from the thermal expansion. Due to the porosity, an overall stress enhance- 
ment occurs in the cathode phase, as can be observed from the overall light 
blue-green color. Figure [TT^c) shows the Von-Mises stress on the YSZ sur- 
face after rotating the volume 180° around the z-axis. Figure [TTj^d) shows 
the Von-Mises stress in the LSC phase. 

Three types of stress enhancements can be noticed from the simulation 
result. At the cathode-electrolyte contacting surfaces, stress is enhanced 
due to the mismatch of thermal expansion and elastic constants between 
the two materials (see the red arrows in Figs. [TT|c) and (d)). The second 
is the concentrated stress observed at the grooves on the electrolyte surface 
(not contacting the cathode) as shown in the white arrows in Fig. [life). 
The third type is the stress concentration effect at the bottlenecks in the 
cathode phase as shown in the green arrows in Fig. [lljd). The simula- 
tion results demonstrate that the smoothed boundary method can properly 
capture the linear elasticity behavior and the geometric effects based on a 
diffuse-interface defined geometry. 

5.4 Phase Transformations in the Presence of a Foreign Sur- 
face 

The Allen-Cahn equation describes the dynamics of a nonconserved order 
parameter, which can be taken as a model for the ordering of magnetic mo- 
ments [26] and diffusionless phase transformations that involve only crys- 
talline order change (26] . It can also be used as a model for evaporation- 
condensation dynamics |26[ [27] . Here, we use the Allen-Cahn equation to 
examine the evaporation of a droplet on a rough surface. The domain pa- 
rameter is given a ripple-like feature as shown in Fig. |i~2| and pre-smoothed 



using Allen-Cahn dynamics, Eq. (58). The droplet phase is placed on top of 



the boundary, and its shape is evolved by the smoothed boundary formula- 



tion of the Allen-Cahn equation, Eq. (49). The simulation is performed in 



two dimensions using parameters Ax = 1, M = 1 and e = 1 with a domain 
size of L x = 100 and L y = 100. 

The evolution of the droplet surface as it evaporates is illustrated in Fig. 
[l2"|a) as a contour ((/> = 0.5) plotted at equal intervals of 270 in dimensionless 
time. The blue to red colors indicate the initial to final stages. As it evolves, 



it is clear that the contact angle is maintained, as shown in Fig. 12 b). The 
dynamics of the motion of the three-phase boundary is interesting in that 
the velocity changes depending on the angle of the surface (with respect to 
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the horizontal axis), which can be inferred from the change in the density 
of the contours. Since the interfacial energy is assumed to be constant, the 
droplet would prefer to have a circular cap shape. However, the contact 
angle imposes another constraint at the three-phase boundary. When the 
orientation of the surface is such that both of these conditions are nearly met, 
the motion of the three-phase boundary is slow while the droplet continues 
to evaporate. When the orientation becomes such that the shape of the 
droplet near the three-phase boundary must be deformed (compared to the 
circular cap), the three-phase boundary moves very quickly. This leads to 
an unsteady motion of the three-phase boundary. On the other hand, at the 
top of the droplet far from the substrate, the curvature is barely affected by 
the angle of the substrate surface; thus, the phase boundary there moves at 
a speed inversely proportional to the radius. 



5.5 Motion of a Droplet due to Unbalanced Surface Tensions 

As another application, we have modeled a self-propelling droplet. Here, two 
different contact-angle boundary conditions are imposed on the right and 
left sides of the droplet placed on a flat surface. The smoothed boundary 



formulation of the Cahn-Hilliard equation, Eq. (53), is used with J n = in 
this simulation. The parameters used are Ax = 1, M = 1, and e = 1, and 
the domain sizes are L x = 240 and L y = 60. The contact angle on the right 
side of the droplet is set to 45 degrees and that on the left side to 60 degrees 
by imposing position dependent boundary conditions. Note that this setup 
is equivalent to the situation in which the wall-environment, droplet-wall 
and droplet-environment surface energies satisfy the conditions of Young's 
equation as 

Iwe — Jwd = Ide cos 60° for the left side, (69a) 



Iwe - Iwd = Ide cos 45° for the right side, (69b) 

where j we , "f w d and 7^ e are the surface energies of the wall-environment, 
droplet-wall and droplet-environment surfaces, respectively, Therefore, this 
model can be used to simulate a case where the surface energies are spatially 
and/or temporally dependent on other fields, such as surface temperature 
or surface composition, as in Ref [64 . 



The evolution of the droplet surface is illustrated in Fig. 13 The droplet 
initially has the shape of a hemisphere, with a 90-degree contact angle with 
the wall surface. The early evolution is marked by the evolution of the 
droplet shape as it relaxes to satisfy the contact-angle boundary condition, as 



26 



seen in Fig. |l3[a). Then the droplet begins to accelerate. Once the contact 
angle reaches the prescribed value, it is maintained as the droplet moves 
toward the right (see Fig. [l~3^b)). In the steady state, the droplet moves at 
constant speed without other effects present. Such motions of droplets have 
been observed and explained as a result of an unbalanced surface tension 
between the head portion (with a dry surface) and tail portion (with a wet 
surface) due to the resulting spatially varying composition and composition- 
dependent surface energy [53] . 

Figure [14] shows the relaxation of an initially hemispherical droplet on 
an irregular substrate surface in a 3D simulation. The contact-angle bound- 
ary condition imposed at the three-phase boundary is 135 degrees. The 
computational box sizes are L x = L y = 120 and L z = 80. As can be seen, 
the droplet changes its shape to satisfy the imposed contact angle, and the 
droplet evolves to a shape for which the total surface energy is minimized. 
The behavior favoring dewetting imposed by the contact angle {9 > 90°) 
is properly reflected in the lifting of the droplet, as shown in Figs. [I4^a)- 
(c) and (d)-(f). During this relaxation process, the three-phase boundary 
shrinks toward the center as the droplet-wall contacting area decreases, as 
shown in Fig. [l4[a)-(c). 



6 Discussions and Conclusions 

In this paper, we have demonstrated a generalized formulation of the smoothed 
boundary method. This method can properly impose Neumann and/or 
Dirichlet boundary conditions on a diffuse interface for solving partial differ- 
ential equations within the region where the domain parameter tp uniformly 
equals 1. The derivation of the method, as well as its implementation, 
is straightforward. It can numerically solve differential equations without 
complicated and time-consuming meshing of the domain of interest since the 
domain boundary is specified by a spatially varying function. Instead, any 
gird system, including a regular Cartesian grid system, can be used with 
this method. 

This smoothed boundary approach is flexible in coupling multiple dif- 
ferential equations. We have demonstrated how this method can couple 
bulk diffusion and surface diffusion into one single equation while the two 



equations serve as the boundary condition for each other in Section 3.3 
In principle, this method can couple multiple differential equations in dif- 
ferent regions that are defined by different domain parameters. For ex- 
ample, the physics within a domain defined by tpi = 1 is governed by a 
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differential equation Hi. The overall phenomenon will be then represented 
by H = s ^ i ipiHi, where the subscript denotes the i-th domain, and 

ipi = 1 represents the entire computational box. When sharing the dif- 
fuse interfaces between domains, the physical quantities can connect to one 
another as boundary conditions for each equation in each domain. There- 
fore, this method could be used to simulate coupled multi-physics and/or 
multiple-domain problems, such as fluid-solid interaction phenomena or dif- 
fusion in hetero-polycrystalline solids. 

We have also demonstrated the capability of applying the smoothed 



boundary method to moving boundary problems in Section 5.2 When the 



locations of domain boundaries are updated by a phase-field type dynam- 
ics such that the domain parameter remains uniformly at 1 and on each 
side of the interface, the smoothed boundary method can be conveniently 
employed to solve differential equations with moving boundaries. 

In addition to Neumann and Dirichlet boundary conditions, we have also 
shown the capability of the smoothed boundary method for specifying con- 
tact angles between the phase boundaries and domain boundaries (Sections 



5.4 and 5.5). This type of boundary condition is difficult to impose using 



conventional sharp interface models. 

Although the smoothed boundary method has many advantages, as 
shown in Section |3.2| the nature of the diffuse interface inevitably intro- 
duces an error proportional to the interfacial thickness since we smear an 
originally zero-thickness boundary into a finite thickness interface. Another 
error results from the resolution of the rapid transition of the domain param- 
eter across the interfacial region. When numerically solving the smoothed- 
boundary formulated equations, properly capturing the gradient of the do- 
main parameter across the interface becomes very important. From our 
experience, at least 4~6 grid points are necessary to resolve the diffuse in- 
terfaces such that the errors are controlled. Moreover, when solving time 
dependent equations, one singularity occurs because of the terms of 1/tp and 
1/ip 2 for imposing Neumann and Dirichlet boundary conditions, respectively. 
In practice, cutoffs at small tp values are necessary to avoid numerical insta- 
bilities. These cutoff values can be smaller as the diffuse interface is better 
resolved, i.e., by using more grid points across the interface. However, only 
a small number of grid points will be used across the interface for compu- 
tational efficiency. In our simulations, when 4~6 grid spacings are used for 
the interfacial regions, the cutoff values are around 1 x 10~ 6 ~1 x 10~ 8 for 
the Neumann boundary condition and 1 x 10~ 2 ~ 1 x 1CP 4 for the Dirichlet 
boundary condition to maintain numerical stability while keeping the er- 
rors reasonably small. On the other hand, when solving time independent 
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equations such as the mechanical equilibrium equation and the steady state 
diffusion equation, there are no singular terms in the equations. The cutoff 
value is simply used to avoid the singularity of the matrix solver. In this 
case, the cutoff value can be as small as the order of numerical precision, 
such as 1 x 1CP 16 . All of these numerical instability and error behaviors re- 
quire more systematic and theoretical studies; thus, the interfacial thickness 
and resolution should be optimized for future works. 

Based on the general nature of the derivation, the smoothed boundary 
method is applicable to generalized boundary conditions (including time- 
dependent boundary values that are important for simulating evolution of 
many physical systems). Since the domain boundaries are not specifically 
defined in the smoothed boundary method, this method can be applied to 
almost any geometry as long as it can be defined by the domain parameter. 
This is very powerful and convenient for solving differential equations in 
complex geometries that are often difficult and time-consuming to mesh. As 
three-dimensional image-based calculations are more prevailing in scientific 
and engineering research fields [65] in which voxelated data from serial scan- 
ning or sectioning are often utilized and are difficult to render as meshes, the 
smoothed boundary method is expected to be widely employed to simulate 
and study physics in complex geometries defined by 2D pixelated and 3D 
voxelated data with a simply process of smoothing the domain boundaries. 
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Figure 1: (a) The conventional sharp interface description of a domain with 
a zero-thickness boundary, (b) The diffuse interface domain and boundary 
defined by a domain parameter, ip. (c) The inward normal vectors defined 
by VV>- 
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Figure 2: Demonstration of the smoothed boundary method on the ID dif- 
fusion equation. The red dashed line is the domain parameter, and the blue 
lines are the concentration profiles taken at different times. The Neumann 
BC is imposed at the right boundary, and the Dirichlet BC is imposed at 
the left boundary. 




Figure 3: The concentration profiles (a) for Df, = 1, k = 2.1, D s = 10 and 
(b) for Df) = 1, k = 1000, D s = 10 obtained by the smoothed boundary 
method (top) and the finite difference method with sharp interface model 
(bottom). The top region with constant blue color is outside of the solid, 
while the solid white lines indicate the solid surface. 
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Figure 4: The concentration along the line at r = 0, r = 2R/3 and r = R 
for (a) D b = 1, k = 2.1 and £> s = 10; and (b) D b = 1, k = 1000 and 
D s = 10. The solid lines are the solutions from the finite difference method 
with sharp interface model. The circular markers are the solutions from 
the smoothed boundary method with 4.5 grid spacings across the interface, 
Ax = 1/60 and R = 60 Ax = 1. The solid square markers are the smoothed 
boundary solutions with 4.5 grid spacings across interface, Ax = 1/30 and 
R = 30 Ax = 1. To clearly illustrate the concentration profile for the low 
surface reaction case, a magnified view is provided in (a). 
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Figure 5: (a) Allen-Cahn type phase transformation with a 60° contact-angle 
BC. (b) Magnified view of the order parameter profile at the three-phase 
boundary, corresponding to t = 2.4 x 10 3 in (a), (c) Cahn-Hilliard type 
phase transformation with a 120° contact-angle BC. (d) Magnified view of 
the order parameter profile at the three-phase boundary, corresponding to 
t = 1.3 X 10 5 in (c). The imposed contact angles can be clearly verified in 
(b) and (d). The field order parameters in the region of < 0.5 have no 
physical significance. For Cahn-Hilliard case, the field order parameter is 
conserved in the region of -0 > 0.5. 
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Figure 6: Steady state concentration for D\, = 1, k = 0.1 and D s = in a 
real cathode complex microstructure: real part (a) and imaginary part (b). 
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Figure 7: The concentration for Df, = 1, n = 2.1 and D s = 10 at intermedi- 
ate time (a) t = 2 x 10~ 2 , and (b) t = 4.9 x 10 _1 in a real cathode complex 
microstructure. Figures (c) and (d) are the magnified views of (a) and (b). 
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Figure 8: (a) Micro-structure of the electrolyte and cathode in SOFC, in 
which cyan and yellow colors indicate the electrolyte (YSZ) and cathode 
(LSC) phases, respectively, (b) The steady-state concentration in the com- 
plex cathode phase for Df, = 1, k = and D s = with fixed C = 1 at the 
irregular LSC- YSZ contact surface. 
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Figure 9: Left column: Normalized concentration profiles (to the lattice site 
density) taken at difference times. Right column: Velocity fields taken at 
different times corresponding to the left column. Black and gray arrows de- 
note the flow inside and outside the material, respectively. The flow outside 
of the material has no physical significance to the shape change. 
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Figure 10: Top row: snapshots of the fast diffuser mole fraction taken at 4 
different times (a)— (d) from initial to final stages. Bottom row: snapshots 
of vacancies mole fraction normalized to the equilibrium value taken at 4 
different times, (e)— (h) corresponding to (a) — (d). The white solid contour 
lines indicate the locations of the rod and void surfaces. 
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Figure 11: (a) Solid phase containing cathode (yellow) and electrolyte (cyan) 
in SOFC. (b) Von-Mises stress in the entire solid phase due to thermal 
expansion, (c) Von-Mises stress in the electrolyte phase, (d) Von-Mises 
stress in the cathode phase. 
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Figure 12: (a) The dynamics of evaporation of a droplet on a rough surface 
(dashed line) governed by Allen-Cahn dynamics. The contact angle between 
the droplet and the surface is imposed at 135 degrees. Solid curves with 
various colors represent the profile of the droplet at different times. The 
outermost blue line represents the initial state, and the innermost represents 
the final state (taken before complete evaporation in the simulation); each 
line is plotted at every time interval of 270 in dimensionless time. The 
velocity of the three-phase boundary is greatly affected by the surface profile, 
(b) A magnified view of the three-phase boundary to show the contact angle 
is accurately set (the angle made by the thin black lines is 135 degrees, (c) 
The order and domain parameters are shown to illustrate the diffuse nature 
of the interface and boundary (taken at t = 270). 
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Figure 13: A self-propelling droplet driven by unbalanced surface tensions. 
The evolution is modeled by the Chan-Hilliard equation with two differ- 
ent contact-angle BCs on each side of the droplet, (a) The droplet shape 
changes during the relaxation period. The color contours are plotted at time 
intervals of 2 x 10 4 in dimensionless times, (b) The droplet motion along 
the substrate surface. The color contours are plotted at time intervals of 
1 x 10 5 in dimensionless time. The droplet moves at constant speed in the 
steady state. 



45 



1 ^[ 




(d) 


in 





Figure 14: A droplet relaxing toward the equilibrium shape. The evolution 
is modeled by the Chan-Hilliard equation with a contact angle of 135° to the 
irregular substrate surface: (a) initial (t = 0), (b) intermediate (t = 3 x 10 3 ), 
and (c) equilibrium state (t = 2.35 x 10 4 ). The three-phase boundaries are 
illustrated by the red color. The side views of the droplet are shown in (d), 
(e) and (f) corresponding to (a), (b) and (c), respectively. 
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